In [1]:
import pylab as pl
%matplotlib inline 

import sys
folder =    '/Users/peppe/work' 
if folder not in sys.path:
    sys.path.append(folder)
import PS4C
from astropy import units as u 
from astropy.coordinates import SkyCoord
from astropy.io import ascii 
import healpy as hp 
In [25]:
from IPython.display import HTML

HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
Out[25]:
In [6]:
from PS4C.PS4Cast  import *
import astropy 
from astropy import units as u,  constants as C

dir_ps='/Users/peppe/work/PS4C/'

PICASSO inpainting SPASS maps


Giuseppe Puglisi , Xiran Bai , Yuuki Omori Sept. 30 ,2019

We implemented the Python Inpainter for Cosmological and AStrophysical SOurces (PICASSO) which is a package aiming at inpaint any kind of maps encoding either Gaussian (e.g. CMB ) or Non-gaussian signals (e.g. Gal. Foregrounds). We made use of several state of the art techniques already adopted for image reconstruction. Most of'em are based on Convolutional Neural Networks (CNN), one of the most promising is the DCGAN . In this posting gonna show the inpainting performances on SPASS TQU maps with DCGAN trained on CMB maps (notice that this isn't logically correct: the training with non-gaussian signal is still ongoing but we oversees to have the trained weights by the end of the week).

How many sources do we expect to detect in the QU SPASS maps?

I run point source forecasts with PS4Cast package (Puglisi + 2018) with the SPASS specs, assuming a detection threshold to be $7 \sigma$ ( this threshold is set by the source detection algorithm ) .

In [7]:
nu=2.3 
sens= .023
fwhm= 8.9 
fsky= 0.3

SPASS = Experiment(ID='SPASS', sensitivity= sens, frequency=nu , fwhm=fwhm , fsky=fsky,nchannels=1,
                units_sensitivity='Jy',units_beam='arcmin')

fc =Forecaster(SPASS,  ps4c_dir=dir_ps, sigmadetection=7 )
fc.forecast_pi2scaling(verbose=False )
fc ()
fc.print_info() 
/Users/peppe/miniconda/envs/pyML/lib/python3.7/site-packages/numpy/ma/core.py:2139: RuntimeWarning: invalid value encountered in greater_equal
  condition = (xf >= v1) & (xf <= v2)
/Users/peppe/miniconda/envs/pyML/lib/python3.7/site-packages/numpy/ma/core.py:2139: RuntimeWarning: invalid value encountered in less_equal
  condition = (xf >= v1) & (xf <= v2)
/Users/peppe/work/PS4C/pointsource_utilities.py:271: RuntimeWarning: invalid value encountered in log
  return np.log(S1/S2)/ np.log(v1/v2)
/Users/peppe/work/PS4C/Stats.py:44: RuntimeWarning: invalid value encountered in log
  return A* 1./(X*sigma *np.sqrt(2*pi)) * np.exp(- (np.log(X/Xm ))**2 /2./sigma**2)
/Users/peppe/work/PS4C/pointsource_utilities.py:514: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and calling the integrator 
  on the subranges.  Perhaps a special-purpose integrator should be used.
  dNdP[i],err= integrate.quad(integrand,a=P, b=S.max().value )
Star Forming galaxies not accounted given the range of frequency you are forecasting. 
==========================================================================================
SPASS	Specifics

Frequency	...... 2.3 GHz
Flux limit	...... 0.023 Jy
Resolution	...... 8.9 arcmin
# channels	...... 1
Fraction of sky	...... 0.3
Beam angle	...... 7.594455960608548e-06 sr
//////////////////////////////////////////////////////////////////////////////////////////
Forecasted quantities for Radio Galaxies 

Frequency 	#sources[S,P]	 Confusion	 <Pi> 	 <Pi^2>x1e3 	D^TT(lensing)	D^BB(lensing)

2.3 GHz	11909	73	43.0716 mJy	4.21	1.89	  9.47232e+09 uK2	  8.93545e+06 uK2
==========================================================================================
Forecasted quantities for Dusty Star Forming  Galaxies 

Frequency 	#sources[S,P]	    <Pi^2>x1e3 	D^TT(lensing)	 D^TT (clustering, lens.)	D^BB(lensing)

==========================================================================================

Notice that the number quoted above are very close to the one reported in Meyers et al. 2017 ( detection of total flux point sources in SPASS survey , no polarized flux). They 've detected about twice the sources we expect, ~23,000 sources since the detection threshold is lower ( $5\sigma$ ). In polarized flux we expect a definitely smaller number of sources ~ 70, with completeness in polarization at $100 $ mJy.

Performing the Source detection algorithm on the SPASS polarization maps

Since we are very interested in polarization, i run the Matched Filter in the polarization amplitude map and if a source is detected with signal-to-noise ratio >7 is included in the detection catalogue.

In [4]:
blobmap= hp.read_map('/Users/peppe/work/heavy_maps/Pfiltered_spasswhole.fits', verbose=False )
spassmap =hp.read_map('/Users/peppe/work/heavy_maps/spass_dr1_1902_healpix_Tb.tqu.fits', field=[0,1,2], verbose=False)

wmask = hp.read_map('/Users/peppe/work/heavy_maps/apodized_spass_footprint.fits.gz', verbose=False) 


nside=hp.get_nside(wmask)
In [5]:
obspixs = pl.ma.masked_not_equal(wmask,0. ).mask
snrmap =blobmap   /blobmap[obspixs].std() 
#print (blobmap[obspixs].std() )
exclude_bright_source = pl.where(snrmap[obspixs]<50)[0]
pl.figure(figsize=(15,15)) 

snrmap =blobmap *wmask /blobmap[obspixs][exclude_bright_source  ].std() 
#print  (( blobmap[obspixs].std() -blobmap[obspixs][exclude_bright_source  ].std()  ) /blobmap[obspixs].std() )
beampix = hp.nside2resol(nside,arcmin=True) 

conv_factor= PS4C.pointsource_utilities.Krj2brightness(freq=2.3*u.GHz, 
                                                     theta_beam= 8.9  *u.arcmin,
                                                     temp= 1. *u.K ) 

hp.gnomview(  pl.sqrt(spassmap[2]**2 +spassmap[1]**2)  
            ,   coord=['G'] ,reso=6
            ,     ysize=300 ,max=0.05,
              xsize=700, rot=[0,-60]  ,sub=311, title='SPASS Polarizat. Amplitude', notext=False, unit='K') ;

obspixs2 = pl.ma.masked_not_equal(wmask,0. ).mask

rms = (blobmap[obspixs2]   *conv_factor  ).std().value
hp.gnomview(  blobmap  *conv_factor 
            ,   coord=['G'] ,reso=6
            ,     ysize=300 ,
              xsize=700, rot=[0,-60],sub=312,  min=-5*rms, max=5*rms, cmap=pl.cm.Greys_r, title='Match Filter map ',
            unit='mJy',notext=False)  

hp.gnomview(  snrmap  
            ,   coord=['G'] ,reso=6
            ,     ysize=300 ,
              xsize=700, rot=[0,-60],sub=313,  min=0, max=5, cmap=pl.cm.Blues, title='SNR map ',notext=False) 
/Users/peppe/miniconda/envs/pyML/lib/python3.7/site-packages/healpy/projaxes.py:1196: MatplotlibDeprecationWarning: 
The iterable function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use np.iterable instead.
  if matplotlib.cbook.iterable(value):
/Users/peppe/miniconda/envs/pyML/lib/python3.7/site-packages/healpy/projaxes.py:1155: MatplotlibDeprecationWarning: 
The iterable function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use np.iterable instead.
  if matplotlib.cbook.iterable(value):
/Users/peppe/miniconda/envs/pyML/lib/python3.7/site-packages/healpy/visufunc.py:552: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  pylab.draw()
/Users/peppe/miniconda/envs/pyML/lib/python3.7/site-packages/IPython/core/pylabtools.py:128: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  fig.canvas.print_figure(bytes_io, **kw)

Constructing the Catalogue

In [6]:
from sklearn.neighbors  import  DistanceMetric
Haversine_distance = DistanceMetric.get_metric('haversine') 
 
def peakfinder(pix, snrmap):
    nside= hp.get_nside(snrmap)
    
    tpeak  = snrmap[pix].max()
    neigh = hp.get_all_neighbours(lonlat=True,nest=False,nside=nside,theta=pix )

    while tpeak < snrmap[neigh].max():
    
        tpeak = snrmap[neigh].max()
        maskpix =pl.ma.masked_equal(snrmap[neigh], tpeak  ).mask 
        pix, tpeak= peakfinder(  neigh[maskpix], snrmap) 
        #print "peak found", tpeak, pix     
     
    return pix,tpeak 


def sizefinder(peak , C, r,delta,snrmap ):

    r2=hp.nside2resol(nside )*delta
    p1 =hp.query_disc( nside, C,  r)
    p2 =hp.query_disc(nside,C,  r2)
    ring =pl.array( list ( set(p2)^ set(p1) ))  
    
    if  snrmap[ring ].mean()<0  :
        #print "size found",pl.rad2deg(r)*60
        return r
    else: 
        delta = delta +1.
        r= sizefinder(peak,C, r2,delta, snrmap  )
        return r 

class Source :
    def __init__(self, lon, lat, size , snr,  pixel ): 
        self.lon =lon 
        self.lat = lat
        coords = SkyCoord( l =self.lon  *u.deg , b=self.lat  *u.deg, frame='galactic')
        self.ra = coords.icrs.ra.deg
        self.dec = coords.icrs.dec.deg
        self.size=size 
        self.snr =snr 
        self.pix=pixel
        self.extended=False 
        
class SourceList():
    def __init__(self) : 
        self.sources=[] 
    
    def add_source (self, source) : 

        if len(self.sources) ==0 : 
            self.sources.append(source)
        else: 
            Lat  =pl.array( [     s.lat[0] for s in self.sources ]  ) 
            Lon  =pl.array( [ s.lon[0] for s in self.sources  ]  ) 
            X=pl.deg2rad( [   (s.lon[0] ,s.lat[0] )  for s in self.sources ]  ) 
            pos = pl.deg2rad( (source.lon, source.lat) ).T
            sigma = source.size   
            dmat = Haversine_distance.pairwise(X, pos) 
            if dmat.min() !=0  : 
                self.sources.append(source)
                
    def savelist (self, filename):
        strings =['GLON' ,'GLAT' ,'RA' , 'DEC' ,'SIZE_ARCMIN', 'SNR', 'EXTENDED', 'MASK_RADIUS' ]
        Lat  =pl.array( [     s.lat[0] for s in self.sources ]  ) 
        Lon  =pl.array( [     s.lon[0] for s in self.sources ]  ) 
        Ra  =pl.array( [     s.ra[0] for s in self.sources ]  ) 
        Dec  =pl.array( [     s.dec[0] for s in self.sources ]  ) 
        size  =pl.array( [     s.size *60 for s in self.sources ]  ) 
        SNR  =pl.array( [     s.snr for s in self.sources ]  ) 
        Ext =pl.array( [     s.extended for s in self.sources ]  ) 
        vals =  [ Lon,Lat, Ra,Dec, size,   SNR, Ext, pl.array(self.maskradii ) ]
        ascii.write(table =  vals,
                         names= strings, output= filename , overwrite=True )
        
In [7]:
SNRthreshold=7

blobs=pl.where(snrmap*wmask  >= SNRthreshold )[0]
nside= hp.get_nside(snrmap)
mask=pl.ones_like(snrmap)
r1=hp.nside2resol(nside ) 
beamfwhm =8.9 

SList=  SourceList() 
pmap =pl.sqrt( spassmap[1]**2 +spassmap[2]**2)

for  ipix  in  (  blobs ):
    peakpix,  snrpeak = peakfinder(ipix,snrmap)
    vec  =   pl.array( hp.pix2vec(ipix=peakpix,nest=False, nside=nside)   )
    lon, lat    =   ( hp.vec2ang( vectors=vec, lonlat=True  ) ) 
    size =pl.rad2deg( sizefinder(peakpix, vec, r1, 2 ,snrmap) ) #degree 
    s=Source(lon,lat,size , snrpeak, peakpix  )
    
    if size >  beamfwhm/60: 
        s.extended=True 
    SList.add_source(s)
SList.maskradii=[]
for s in SList.sources: 
    vec  =   pl.array( hp.pix2vec(ipix=s.pix,nest=False, nside=nside)   )
    hole = hp.query_disc(nside=nside, vec=vec, radius=pl.deg2rad( s.size )) 
    T_integrated= pmap[hole].sum() *u.K
    flux =PS4C.pointsource_utilities.Krj2brightness(freq=2.3*u.GHz, 
                                                     theta_beam= s.size  *u.deg,
                                                     temp=  T_integrated).to(u.Jy)
    f = 1. + pl.log10(flux.value)
    if s.extended: 
        radius = pl.deg2rad(s.size *f     )
    else: 
        radius = pl.deg2rad(beamfwhm/60.  )
    SList.maskradii.append(radius)
    disc = hp.query_disc(nside=nside, vec=vec, radius=radius ) 
    mask[disc ] = 0.
        

print ("Found %d point sources w/ SNR> %d" %( len( SList.sources ), SNRthreshold) )
Ext =pl.array( [     s.extended for s in SList.sources ]  ) 
print ("%d are extended  " %( len( Ext[Ext] ) ))
SList.savelist('/Users/peppe/work/inpainting/FG_inpainting/SPASS_pol_pointsources_SNR{}.dat'.format(SNRthreshold))
hp.write_map("/Users/peppe/work/heavy_maps/SPASS_pointsource_mask_matchfilter_SNR{}.fits".format(SNRthreshold), mask, overwrite=True  )
Found 558 point sources w/ SNR> 7
340 are extended  

As you can notice there are ~500 sources most of them are in the Galactic Plane ( in PS4Cast these aren't included ).Thus, we have to apply a mask which selects the Extra-galactic sources (setting a threshold in Galactic latitude $|b|>20 ^{\circ}$, excluding the extended sources ( i.e. the ones whose size is larger than the beam) .

In [26]:
strip = hp.query_strip(nside, pl.deg2rad(70),pl.deg2rad(110) )
strippedmask=wmask*1.
strippedmask[strip]=0. 
ones = np.where( strippedmask>0)[0]
In [27]:
hp.mollview(mask*strippedmask)
hp.graticule()
/Users/peppe/miniconda/envs/pyML/lib/python3.7/site-packages/healpy/projaxes.py:1196: MatplotlibDeprecationWarning: 
The iterable function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use np.iterable instead.
  if matplotlib.cbook.iterable(value):
/Users/peppe/miniconda/envs/pyML/lib/python3.7/site-packages/healpy/projaxes.py:1155: MatplotlibDeprecationWarning: 
The iterable function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use np.iterable instead.
  if matplotlib.cbook.iterable(value):
/Users/peppe/miniconda/envs/pyML/lib/python3.7/site-packages/healpy/visufunc.py:296: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  pylab.draw()
/Users/peppe/miniconda/envs/pyML/lib/python3.7/site-packages/healpy/visufunc.py:1339: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  pylab.draw()
/Users/peppe/miniconda/envs/pyML/lib/python3.7/site-packages/IPython/core/pylabtools.py:128: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  fig.canvas.print_figure(bytes_io, **kw)
0.0 180.0 -180.0 180.0
The interval between parallels is 30 deg -0.00'.
The interval between meridians is 30 deg -0.00'.

Above the point source mask excluding the low Gal. latitude sources.

In [28]:
tab = ascii.read('/Users/peppe/work/inpainting/FG_inpainting/SPASS_pol_pointsources_SNR7.dat')
In [29]:
Nentries =tab['RA' ].shape[0]
extended =  tab['SIZE_ARCMIN'] >=10.

Selecting unresolved sources at $|b|> 20 ^{\circ}$

In [30]:
s =0 
maps=pl.sqrt( spassmap[2]**2 + spassmap[1]**2)*( wmask )
masklat= pl.ma.masked_greater(abs( tab['GLAT'] ), 20 )
masksize=pl.ma.masked_less_equal(tab['SIZE_ARCMIN'],11 )
maskand = pl.logical_and ( masklat.mask,  masksize.mask )
print(tab['GLAT'][maskand].shape) 

        
        
for lon ,lat in zip(tab['GLON'][maskand]  , tab['GLAT' ][maskand]) :
    pl.figure()
    
    hp.gnomview(maps , rot=[lon, lat], 
                   xsize=128, reso=3.5,sub=121, min=0 )
    hp.gnomview(( snrmap ), rot=[lon,lat],min=0., max=7, xsize=128,
                cmap=pl.cm.Blues, reso=3.5,sub=122) 
(45,)
/Users/peppe/miniconda/envs/pyML/lib/python3.7/site-packages/healpy/visufunc.py:552: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  pylab.draw()
/Users/peppe/miniconda/envs/pyML/lib/python3.7/site-packages/ipykernel_launcher.py:11: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  # This is added back by InteractiveShellApp.init_path()

Notice that we end up considering 45 sources, which is very close to our expectations.

Inpainted Maps

The masked area is centered in the detection coordinates and a masking radius of 30 arcminutes (3x SPASS beam).

Now let's see what Picasso can do, let's inpaint'em all!

In [4]:
ra,dec= pl.loadtxt(folder+'/inpainting/FG_inpainting/SPASS_sources_snr7_glon_glat.dat')
for i in range(ra.shape[0]):
    pl.figure()
    im = pl.load(folder+'/inpainting/FG_inpainting/stacks/SPASS/singlestacks/Q_{:.5f}_{:.5f}_masked.npy'.format(ra[i],dec[i] )) 
    inp = pl.load(folder+'/inpainting/FG_inpainting/outputs/Contextual-Attention/Q_{:.5f}_{:.5f}.npy'.format(ra[i],dec[i] )) 
    pl.subplot(121)
    pl.title('Masked SPASS')
    pl.imshow(im);pl.colorbar(orientation='horizontal')    
    pl.subplot(122)
    pl.title('Inpainted SPASS')
    
    pl.imshow(inp);pl.colorbar(orientation='horizontal')
    
    
    
/Users/peppe/miniconda/envs/pyML/lib/python3.7/site-packages/ipykernel_launcher.py:3: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  This is separate from the ipykernel package so we can avoid doing imports until

These are 6x6 deg flat sky Q maps (128 x 128 pixels ) extracted from the SPASS map and centered at the source location (sizes are similar to the maps shown above). As already stated above, this can be definitely improved once we 've the trained weights on Synchrotron simulated maps. Overall looks very promising!

In [ ]: